Convective Heat Transfer in the Reusable Solid Rocket Motor of the Space 

Transportation System 

Rashid A. Ahmad 

Motor Performance Department, Gas Dynamics Section, M/S 252 
Science and Engineering 

ATK Thiokol Propulsion, A Division of ATK Aerospace Company Inc., Brigham City, Utah 

84302 

This simulation involved a two-dimensional axisymmetric model of a full motor initial grain 
of the Reusable Solid Rocket Motor (RSRM) of the Space Transportation System (STS). It was 
conducted with CFD commercial code FLUENT.® This analysis was performed to. a) maintain 
continuity with most related previous analyses, b) serve as a non-vectored baseline for any three- 
dimensional vectored nozzles, c) provide a relatively simple application and checkout for various 
CFD solution schemes, grid sensitivity studies, turbulence modeling and heat transfer, and d) 
calculate nozzle convective heat transfer coefficients. The accuracy of the present results and the 
selection of the numerical schemes and turbulence models were based on matching the rocket 
ballistic predictions of mass flow rate, head end pressure, vacuum thrust and specific impulse, 
and measured chamber pressure drop. Matching these ballistic predictions was found to be good. 
This study was limited to convective heat transfer and the results compared favorably with 
existing theory. On the other hand, qualitative comparison with backed-out data of the ratio of 
the convective heat transfer coefficient to the specific heat at constant pressure was made in a 
relative manner. This backed-out data was devised to match nozzle erosion that was a result of 
heat transfer (convective, radiative and conductive), chemical (transpirating), and mechanical 
(shear and particle impingement forces) effects combined. 


Presented as Paper AIAA-200 1-3585, AIAA/ASME/SAE/ASEE 37th Joint Propulsion 
Conference, Salt Lake City, UT, July 8-11, 2001. 

*Sr. Principal Engineer, Associate Fellow AIAA. 

®2002 ATK Thiokol Propulsion, a Division of ATK Aerospace Company Inc., Published by the 
American Institute of Aeronautics and Astronautics 


l 


Nomenclature 


o 2 

A area, m" (ft ) 

a empirical constant (Eq. 3) 

Cf friction coefficient 

C p , C v specific heats at constant pressure and volume, respectively, J/kg-K (Btu/lb m -R) 

COD coefficient of determination of a correlation 
d diameter, m (in.) 

2 2 

g acceleration due to gravity, m/s (ft/s ) 
h convective heat transfer coefficient, W/irf-K (Btu/hr-ft'-R) 

I turbulence intensity, w’/w«, (%) 

K an acceleration parameter 

k gas thermal conductivity, W/m-K (Btu/hr-ft-R) 

/ mixing length constant usually known as VonKarman length scale, m (ft) 
m mass flow rate, kg/s (lbm/s) 

MW molecular weight, kg/kgmole (lb/lbmole) 

M Mach number 

n burning rate pressure exponent 

Po , P total and static pressure, respectively, Pa (psia) 

Pr Prandtl number 

Pr, turbulent Prandtl number 

q wall heat flux, W/m 2 (Btu/hr-ft 2 ) 

Re(x) local axial Reynolds number based on distance along the nozzle wall, pjx) «„ (x) x w / 
JJbdx) 
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Red(x) local axial Reynolds number based on diameter of the nozzle, pjx) u„ (x) d(x) / pjx) 
recovery factor 

f propellant bum rate, m/s (in./s) 

r radius, m (in.) 

S source term (Eq. 1) 

To, T total and static temperature, respectively, K (R) 
f + non-dimensional temperature in wall coordinates 

и, v axial and radial velocity components, m/s (ft/s) 

m<„ velocity magnitude of u and v at the motor centerline, m/s, ft/s 

m + non-dimensional velocity in wall coordinates 

m ' velocity fluctuation, m/s (ft/s) 

m t shear velocity, m/s (ft/s) 

w propellant weight flow rate, N/s (lbf/s) 

x, r axial and radial axes, m (in.) 

(x), (r) function of axial and radial directions, respectively 

y + non-dimensional distance from the wall in wall units 

a augmentation/de-augmentation factor in propellant mass flux 

y ratio of specific heats, Cp/C v 

к, e turbulence kinetic energy (m 2 /s 2 (ft 2 /s 2 )) and its rate of dissipation (m7s7s (ft7s7s)), 

respectively 

p, v dynamic (N-s/m 2 (lb m /ft-s)) and kinematic (m7s (ft7s)) viscosity, respectively 
p gas density, kg/m 3 (lbm/ft 3 ) 

7 7 

Xw wall shear stress, N/m" (lbf/ft') 
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Subscripts 


aw adiabatic wall or recovery temperature 

e at nozzle exit 

/ face 

g gas 

h hydraulic 

o chamber conditions 
p propellant 

ref reference 

s static 

t turbulent 

vac vacuum 

w wall 

oo conditions at motor centerline 

Superscripts 
* throat conditions 

” flux 


mean and time averaged 



Introduction 


A computational fluid dynamics (CFD) simulation involved a two-dimensional axisymmetric 
model of a full motor initial grain of the Resusable Solid Rocket Motor (RSRM) was conducted 
with CFD commercial code FLUENT. 1, 2 This analysis was performed to a) Maintain continuity 
with the most related previous analyses 3-7 of this motor, b) Serve as a non-vectored baseline for 
any three-dimensional vectored nozzles; c) Provide a relatively simple application and checkout 
for various CFD solution schemes, grid sensitivity studies, turbulence modeling/heat transfer, 
etc.; and d) Calculate the nozzle convective heat transfer coefficients. Nozzle heat transfer is of 
interest because the convective heat transfer coefficient to the specific heat at constant pressure 
ratio ( h/Cp ) is usually used as input in the Charring and Material Ablation (CMA) code 8 for 
nozzle erosion predictions. 

The following are the four most related studies of this motor. First, Golafshani and Loh 3 
conducted a time-dependent, axisymmetric numerical solution of the Navier-Stokes to analyze 
the viscous coupled gas-particle non-reacting flow in solid rocket motors. The solution assumed 
laminar internal flow. Second, Loh and Chwalowski 4 used particles of diameters of 1 to 100 |xm 
in converging-diverging nozzles and a mass loading of 28.8%. Acceleration between lg and 3g 
had a minimal effect on the particles’ behavior in the nozzle. Third, Whitesides et al. 5 conducted 
a two-dimensional axisymmetric two-phase flow analysis in the RSRM. The overall objective 
was to determine the structure of the flow field in the recirculation region underneath the 
submerged nozzle nose and to define the gas flow and particle impingement environments along 
the surface of the aft case dome insulation. It was concluded that particles were impacting the 
area underneath the nozzle nose and forming a sheet of molten aluminum oxide or slag. The 
sheet flows afterwards, along the underneath nozzle nose surface as is the direction of the near 
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surface velocity vector during the last half of motor bum. This slag layer is then sheared from 
the nozzle cowl/boot ring surface and impacts the aft dome case insulation at the location of 
severe erosion. Fourth, Laubacher 6 conducted two-dimensional axisymmetric analysis to 
compute chamber pressure drop in the RSRM. The walls were assumed to be adiabatic and 
utilized the standard k-£ turbulence model and a coupled solver using in-house code, SHARP. 
These studies were conducted for the RSRM and no attempt was made to calculate the 
convective heat transfer coefficients. The fifth study 7 has the details of this study. It involved 2- 
D axisymmetric and 3-D vectored nozzles. Furthermore, it involved two-phase flow where slag 
concentration, accretion rates, and particle trajectories were calculated. At this time, it is not 
feasible how to augment convective heat transfer by effects of particle impingement and 
trajectories. 

Related convective heat transfer studies are given in Refs. (9-16). When considering 
convective heat transfer in solid rocket motors, surface temperatures and heat fluxes are high and 
very difficult to measure. Ablative materials are used to dissipate and inhibit heat transfer by 
erosion and transpiration. For the lack of reliable thermal conditions, the nozzle wall was usually 
assumed to be adiabatic in CFD calculations. On the other hand, CFD calculations (velocity, 
density, pressure, temperature, viscosity, etc.) and geometry enable someone to calculate heat 
transfer. 9 It is usually estimated 9 using three well-known methods. They are the modified 
Reynolds’ analogy 17 for laminar flow over a flat plate, Dittus-Boelter correlation 17 for fully- 
developed turbulent pipe flow, and the Bartz 10 correlation for nozzle flows. Bartz 10 extended the 
well-known Dittus-Boelter correlation 17 for turbulent pipe flow to account for mass flux and 
variations in velocity and temperature. Back et al. 11, 13 15 conducted analytical and experimental 
convective heat transfer studies in the Jet Propulsion Laboratory (JPL) nozzle. Moretti and 
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Kays 12 conducted experimental convective heat transfer to an essentially constant property 
turbulent boundary layer in a two-dimensional channel for various rates of free-stream 
acceleration. Back et al. 14-15 and Moretti and Kays 12 found that acceleration causes a depression 
in heat transfer rate below what would be predicted assuming a boundary-layer structure such as 
obtained for constant free-stream velocity. They attributed it to re-laminarization of the turbulent 
boundary layer. They further state, it is by no means obvious that the same acceleration 
parameter applicable to an axisymmetric flow. Wang 16 focused on the capability of general- 
purpose CFD codes in predicting convective heat transfer coefficients between a fluid and a solid 
surface. Effects of various parameters such as grid resolution, turbulence models and near-wall 
treatments, as well as numerical schemes on the accuracy of predicted convective heat transfer 
were studied. Test cases included flat plate, pipe flow, JPL nozzle, and impinging jets. 

The attributes of this two-dimensional (2D) axisymmetric analysis are: 

• Significant effort was made to assess grid sensitivity and grid consistency with turbulence 
models. 

• Verifying flow/thermal solution quality represented by vs. u + vs. y + and t + vs. y + against the 
usual velocity and thermal laws of the wall for incompressible turbulent flows, respectively. 

• Calculations of nozzle convective heat transfer, including the assessment of turbulent 
boundary layer re-laminarization. 

Discussion of Modeling Approach 

Governing equations, geometry parameters, operating and predicted ballistic conditions, gas 
thermophysical properties, grid density and turbulence modeling, boundary conditions, 
computational schemes and numerical convergence (residuals) are discussed next. 

Governing Equations: The numerical studies considered the solution of the Navier-Stokes 
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equations, energy equation, the turbulence kinetic energy with its rate of dissipation equations, 
mass transfer and the necessary constitutive equations (ideal gas law, power law for gas thermal 
conductivity and viscosity, etc.). The general governing equation was 

V .{pv <t, -T t V 4>) = S, ( 1 ) 

and the mass conservation equation 

V . {p V) = 0 (2) 

where 0 can be velocity components («, v, w), enthalpy (i) and turbulence quantities (k, e); /’is 
an exchange coefficient for <p\ S# is a source term for 0 per unit volume. 

Geometry Parameters: Table 1 gives a summary of geometry parameters. They are the 

normalized chamber and exit radii along with chamber and exit area ratios. 

Operating and Predicted Ballistic Conditions: Table 1 gives a summary of the ballistic 

prediction parameters for the RSRM that include head end pressure and FSM-9 measured 

chamber pressure drop. 19 The accuracy and the selection of the schemes, models and results are 

based on matching the above ballistic prediction parameters in addition to mass flow rate, and 

1 8 

vacuum thrust and specific impulse. 

Gas Thermophysical Properties: The total head-end pressure given in Table 1 along with 
propellant formulation was used as input to the NASA-Lewis program' or the ODE module of 
the SPP code 21 to obtain chamber gas temperature ( T „ ), dynamic viscosity specific heat at 
constant pressure ( C p ), thermal conductivity ( k u ), and molecular weight. They are given in Table 
1. The local gas dynamic viscosity and thermal conductivity were calculated as a function of 
local temperature as given at the bottom of Table 1. 

Grid Density and Turbulence Modeling: Table 2 gives the turbulence models used with the 
desired values for wall y + and the pertinent results. A coarse and fine grids with quadrilateral 
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cells used in this simulation and are summarized in Table 2. They were designed, solved, and 
iterated on to give the desired values for y + so that consistency with the turbulence models was 
achieved as given in Table 2. All the grids were generated by using GRIDGEN " and made 
orthogonal and smoothed from one domain to another. 

Boundary Conditions: Boundary conditions are applied at the propellant surface, nozzle exit, 
and walls and discussed as follows: 

m M th£. propellant vm/flcg: Mass flux was calculated as a function of the local static pressure 
as 

\v> i \T (3 a ) 

m = a p p a[P s {x,r)\ 


where 



(3b) 


In addition, uniform chamber temperature, flow direction that was normal to the propellant 
surface, an assumed turbulence intensity, 9 (7), and hydraulic diameter ( dh ) were specified. The 
augmentation factor, a, was used as 1 for the propellant except in the head end fin region, where 
it was increased to 4.528 to account for the three-dimensional fins modeled in two-dimensional 
axi-symmetric analysis. 

( 2) At exit: A supersonic boundary condition was utilized where the quantities ( P , T, u, v, k, e) 
were calculated from cells upstream of the exit. The exit pressure, temperature, turbulence 
intensity, and exit hydraulic diameter were specified to start the calculation. The exit pressure 
and temperature were updated as the solution proceeded. 


73) At wall: Three wall boundary conditions are used and discussed as follows: 
(a) Velocity wall boundary condition was assumed to be no slip condition. 
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(b) Thermal wall boundary conditions were assumed for the submerged and the converging- 
diverging part of the nozzle walls. The submerged wall was assumed to be isothermal at 2938.5 
K (5289.3 R). On the other hand, the nozzle wall was assumed to be non-isothermal. The 
surface temperature was taken from Ref. 23 and curve-fitted by this author using 
TableCurve2D 24 as 

a+cx+ex 2 + g x* +ix* +kx 5 (4) 

1 + b x + d x 2 + f x 3 + h x 4 + j x 5 

where x was taken along the nozzle surface and where the coefficients are given as follows: a = 
2789.03, b = - 4.61, c = -13307.21, d= 11.96, e = 32905.45,/= - 6.33, g = - 20712.33, h = - 
0.916, y = 1 174.06 and k = 815.51. The correlation coefficient was calculated to be r - 0.997 . A 
user defined function (UDF) was used to compile the specified surface temperature profile. 
Assuming this surface temperature profile enables the calculation of the heat flux that in tum 
enables the calculation of the convective heat transfer coefficient depending on the assumption of 
a reference temperature. Calculated centerline and recovery temperatures are also shown and 
will be discussed shortly. 

(c) Near wall turbulence treatment: Two equation turbulence models were used. They 
are the standard k-£ and RNG k-£. Near wall treatment involved the standard wall functions and 
two-layer zonal models as described in Table 2. The standard wall functions are given in terms 
of the non-dimensional velocity in wall coordinates (w + ), non-dimensional distance from wall in 
wall coordinates (y ) and the frictional velocity (u T ) are defined' as 



for velocity and 
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(5b) 


t + = 


T w ~T(r) 




(P» U r C p) 


A exit plane 


for temperature at the exit plane. 

In the two-layer zonal model, ' " the wall functions are not used. Rather, the turbulent 
boundary layer is divided into two layers distinguished by a wall-distance-based turbulent 
Reynolds number (Re, = p ic m y / //). The first layer is adjacent to the wall and is called a 
viscous sublayer (viscosity (v) is much larger than eddy diffusivity for momentum (em)) where 
the low standard k-e turbulence model (Low-Re) is used. Only the k equation is solved in the 
viscosity affected region while e is computed from a length scale correlation. The second layer is 
fully turbulent (v « em) and where the high standard k-e turbulence model (High-/?e) is used. 
This model requires finer mesh resolution and therefore larger CPU and memory resources. 
Computational Schemes: The segregated solver in the commercial code Fluent was used. 
Differencing schemes utilized were l sl and 2 nd order Upwind, Power law, and Quick schemes. 
The 1 st Upwind scheme was used to start the problem and the higher order schemes to obtain the 
final results. The 2 nd order Upwind and Quick schemes were found to give similar results in 
terms of mass flow rate and mass imbalance, head-end pressure, chamber pressure drop, and 
maximum Mach number at the nozzle exit. 

Numerical Convergence (Residuals): Numerical convergence was achieved by satisfying four 
requirements in the following sequence. First, the residual error diminished as the number of 
iterations was increased. Second, the profiles of the variables ceased to change, at least 
qualitatively. Third, a first monitor on the total pressure at the propellant surface until the 
average total pressure ceased to change. Fourth, a second monitor on the mass imbalance 


li 


between the inlet (propellant surface) and the outlet (nozzle exit) mass flow rates until it reached 
a small value (10 3 - 10 5 kg/s). 

Results 

Motor chamber pressure drop, turbulence results, convective heat transfer and acceleration 
parameter discussed next. 

Motor Chamber Pressure Drop: References 21, 26 and 27 give the chamber pressure as 6.30 
MPa (913.85 psia), 6.24 MPa (905 psia), and 6.27 MPa (910 psia), respectively. The' total 
pressure used in this study was 6.28 MPa (910.78 psia). 

Figure 1 shows the geometry considered and the static pressure distribution in the whole 
motor at 1 s post ignition. Figure 2 shows the submerged cavity and its location at prior to 
ignition relative to the nozzle. Figure 3 shows the local axial static pressure along the centerline 
of the RSRM chamber. The axial coordinate, x, started at the head end and ended at the nozzle 
entrance and was measured along the centerline. The interest was to match the calculated motor 
chamber pressure drop against measured data from static tests of QM-7 and QM-8 (Ref. 27) and 
FSM-9 (Ref. 19). For the cases with the default value for C M (Model 1 and 2a as given in Table 
2), the viscosity ratio (effective (laminar plus turbulent) to laminar) was limited to 10,000 to 
match the total pressure of 6.28 MPa (910.78 psia). The chamber pressure drop was calculated to 
be 1.38 MPa (200 psia). In Model 2b of Table 2, the viscosity ratio was not limited, but C M was 
reduced from the default value of 0.09 to 0.055." Variable C M was suggested by modelers and is 
well substantiated by experimental evidence." For example, C M is found to be around 0.09 in the 
inertial sublayer of equilibrium boundary layers, and 0.05 in a strong homogeneous shear flow. 
This was applied in Model 2b of Table 2. The motor chamber pressure drop was calculated to be 
1.23 MPa (178 psia). A better agreement in chamber pressure drop (1.17 MPa (170 psia)) was 
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achieved on a finer grid using Model 2c of Table 2. On the other hand, the calculated values 7 of 
y + were in the range between 0.5 and 8 and with some scattering. The scattering is because of 
the high values of aspect ratio of cells. In addition, it is suggested 1 ' 2 to locate the first grid point 
adjacent to the wall at y + = 1. Therefore, the results in this study are reported as calculated on a 
coarse grid with y + consistent with the standard wall functions. The results of all turbulence 
models used will be given when compared. Otherwise, the results of Model 2b of Table 2 will be 
given. 

Also shown are the measurements of the local static pressure along QM-7 and QM-8 test 
motors. In addition, chamber pressure drop was measured to be 1.14 MPa (165 psia) in FSM-9 
static test. 19 These measurements were made between the head end and in the vicinity of joint 2 
in the submerged region. The above agreement is acceptable. 

Turbulence Quantities: Two turbulence models were utilized in this study and given in Table 
2b. The calculated values of y + on the fine grid 7 were in the range between 0.5 and 8 and with 
some scattering. The scattering is because of the high values of aspect ratio of cells and the 
desired wall y + were not fully achieved. The size of this motor is too large and axial refinements 
were not pursued because they are cost and time prohibitive. The height of the first cell adjacent 
to the nozzle wall can be estimated as: 


A = 2 (r.-r) = 


2v . 

— y 

U r 


(6) 


and shown in Fig. 4. To obtain y + < 1 would require grid refinement in the radial direction. 
Additional refinement in the axial direction is required to keep the cell aspect ratio acceptable 
numerically and to avoid cells of negative areas. Thus large number of cells (in millions) are 
required for this large motor. It was found that large cell aspect ratio would generate 
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questionable scattering in wall y + and thus heat transfer coefficient. On the other hand, 
calculated chamber pressure drop calculated on the finer grid compared more favorably than the 
coarse grid. Since the interest here is the calculation of the convective heat transfer, the results 
are mainly given on the coarse grid unless stated otherwise. 

Figure 5 shows the calculated wall y + along the converging-diverging part of the nozzle. For 
the nozzle wall values (y + , h, etc.) the axial coordinate, x w , started at the beginning of the 
convergent section of the nozzle and ended at the nozzle exit and was measured along the nozzle 
wall. All the models used give similar wall y + profiles and values as shown in Fig. 5. Since the 
values were between 25 and 130, the first and second turbulence models used were consistent 
with grid resolution. 

Calculated flow quality of the present results is shown in terms of y + vs. u and y + vs. r + and 
in comparison with the velocity and thermal laws of the wall for incompressible turbulent flow. 
Velocity Law of the Wall: Figure 6 shows the present results at the nozzle exit in comparison 

with the famous incompressible turbulent velocity law of the wall for an external boundary layer, 

28 

i.e. the Spalding’s profile: 


y + =u + + exp(- K B ) 


( + \ , + (jcu + ) 2 {kuJ 

exp(/f u j-1 -ku - - 


(7a) 


where k and B were taken 28 as 0.4 and 5.5, respectively. The main assumptions made in the 
above law were incompressible, negligible stream-wise advection, no axial pressure gradient, and 
no transpiration. 

The first flow cell from the present study was located at y + of 29. Very good agreement with 
Spalding’s profile 28 was achieved around y + = 25. Furthermore, at y + = 100, the incompressible 
velocity law of the wall yielded 16.8 for u + . The present calculations for compressible turbulent 
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flow using Model 1, Model 2a, and Model 2b (Table 2) yielded 16.1, 16.3, and 15.8, 
respectively. They are depressed by -4.17%, -2.98%, and -5.95%, respectively. 

The corresponding compressible turbulent velocity law of the wall is a lot more involved and 
is unwarranted. This above comparison is partially in agreement with conclusions made by 
White 28 regarding the compressible turbulent velocity law of the wall. His conclusions were: 

1) The effect of increasing high Mach number (compressibility, y= 9tu T ) / (2 C p T w )), for a 
given Tw/Taw depress h + below the incompressible turbulent velocity law of the. wall. 
This is the case in Model 1 and Model 2a, but partially with Model 2b. 

2) On the other hand, a cold wall (heat flux, (5=q w w / T w k w u T ) tends to raise w + above the 
incompressible turbulent velocity law of the wall. This is the case partially in Model 2b. 

Temperature Law of the Wall: On the other hand, neither compressible turbulent thermal law 
of the wall exist nor conclusions were made by White." To our dismay, again we need to 
compare with incompressible turbulent flow. Similarly, Fig. 7 shows the present results at the 
nozzle exit in comparison with the temperature law of the wall. The incompressible turbulent 
temperature law of the wall for an external boundary layer was taken as 

t + = Pr y + ; y + <13.2 (7b) 

Pr v + ^ 

t + = 13.2 Pr +— £- In ; y + >13.2 

t 1^13.2 J 

where Pr, and / were taken as 0.85 and 0.41, respectively. Once again, the main assumptions 
made in the above law were incompressible, negligible stream-wise advection, no axial pressure 
gradient, and no transpiration. 

Only a qualitative agreement with the thermal law of the wall is achieved. Similarly, at y + = 
100, the incompressible thermal law of the wall yielded 10.9 for t + . The present calculations 
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using Model 1, Model 2a, and Model 2b (Table 2) yielded 8.5, 8.1, and 6, respectively. They are 
depressed by -22.02%, -25.69%, and -44.95%, respectively. The data shown in Fig. 7 from the 
present study is correlated by the present author to give the following: 

r = [o. 1525 + 1.618 ln(y + )]; 25 < y* < 1000, COD = 0.80 (7c) 


from Model 1 , and 

t* = [- 0.4793 + 1 .709 ln(y + )] ; 25 <y*< 1000, COD = 0.83 (7d) 


from Model 2a, and 

t + = [3.965 + 0.4292 ln(y + )] ; 29 < / < 1000, COD = 0.63 (7e) 


from Model 2b of Table 2. The “COD" is the coefficient of determination of a correlation. The 
higher COD, the higher the quality of the correlation. These correlations apply to compressible 
turbulent flow. 

Convective Heat Transfer: The ratio of the convective heat transfer coefficient to the specific 
heat at constant pressure ( h/C P ) was calculated by eight different methods. The eight methods 
will be given, shown and finally discussed. The first two methods used finite volume CFD 
(Methods 1 and 2), followed by four approximate methods (Methods 3, 4, 5 and 6), followed by 
an integral method (Method 7) and then using backed-out data from measurements (Method 8). 

In Method 1, it was calculated internally using the calculated heat flux based on the difference 
between the local normalized specified wall temperature given by Eq. (4) and shown in Fig. 8 
and the chamber temperature used as a reference temperature, i.e.. 


U) _ q „ U) 
c r ~c„ [T.-TSx)} 


(8a) 


and is shown in Fig. 9. 
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In Method 2, it was calculated using the recovery temperature and defined as 


(9a) 

Mf) _ <?H (*) 

C P C p [T aw (x)-T w (x)] 

and shown in Fig. 9. The recovery temperature 25 is given as 

r„. , (x) s T_ (x) + 95, M [r„ - r (*)] (9b) 


and where i(x) = [Pr,(;c)] l/3 was used for turbulent flow and was calculated locally. The 
temperature at the edge of the boundary conditions has been replaced by the local axial static 
temperature along the motor centerline and shown in Fig. 8. The recovery temperature (T a w. i) 
was calculated at the motor centerline and shown in Fig. 8 to be less than the chamber 
temperature. The chamber ( T 0 ) and recovery temperature (Taw, 2 ) were matched only by taking Pr 
as 1 which yields 1 for SH 2 which corresponds to an ideal situation. Also shown the gas 
temperature ( T g , w ) in the cell adjacent and along the nozzle wall. Good agreement between 
Methods 1 and 2 is shown. 

The four approximate methods (Methods 3, 4, 5 and 6) are discussed next. The third method 
used turbulent flow over a flat plate 17 which is 


= 0.0296 [Re(x)]° 8 Pr(jc) ,/3 


~k(x)~ 

( 1 ) 

X 

l c J 


( 10 ) 


and is not shown in Fig. 9, because it overestimates the heat transfer and does not have the 
correct profile. It appears approximating this large nozzle as a flat plate is not plausible. The 
maximum occurred at the start of the nozzle and decreased along the wall as expected for flow 
over a flat plate. 

The fourth method used the modified Reynolds' analogy for laminar flow over a flat plate 
which is, 
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17 


(11a) 


Re(x) Pr(x) l/3 



{ 1 ) 

X 

l c J 


where the local skin friction is calculated as 


C f (x) = 







(Hb) 


and is not shown in Fig. 9. Again, this method overestimates the heat transfer but has the correct 
profile. In addition, Methods 3 and 4 are not shown to reduce clutter. 

The fifth method used the Dittus-Boelter 17 correlation for fully-developed turbulent pipe flow 


as 


■ 0.023 [Re rf (;t)] 0 8 PrC*)" 


'*(*)" 

r_n 

_d(x)_ 

K Cp j 


( 12 ) 


and is shown in Fig. 9. The exponent n for Pr was taken as 0.3 for cooling ( T w < TJ> and 0.4 for 
heating ( T w > TJ) per Fig. 8. 


The sixth method used the Bartz 10 correlation 


h 6 (x) _ 


C. 


0.026 

(cT 


( ,, 0.2 \ 


pr- c 


Pr 


06 


Jo 


fp..*r(d' r 


V c J 


K r < ) 


f A . ^ 


A(x) 


(T(x) 


( k(x)^ f ' ^ 


^d(x) j 


v c v 


(13a) 


and shown in Fig. 9. It was based on a similar correlation of fully developed turbulent pipe flow 
heat transfer, as given previously. The terms in the bracket are constant for a single nozzle. The 
subscript “o”, on the second term in the bracket, signifies properties are to be evaluated at 
chamber conditions. Prandtl number at chamber conditions was evaluated 10 using 

4 y (13b) 


Pr = 


9 ;y-5 


and was calculated to be 0.843. The ratio of specific heats is calculated from C P and the gas 
constant. Using lower values for Pr (= 0.5) as in this study, would overestimate h/C p as 
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calculated by Eq. (13a). The characteristic exhaust velocity was calculated as 


{yRTj 2 



(13c) 


The first term outside the bracket in Eq. (13a) is a function of nozzle local area ( A(x ) = (tc/4) 
x )). The second term outside the bracket, o(x), is a dimensionless factor accounting for 
variations of density and viscosity across the boundary layer. It is shown in Bartz 10 to be 


a(x) = - 


(13d) 


I JM [, + Zri« 


ll 

0.8 — — 
5 

1 

1 <N 

+ 





where co is taken to be 0.6. This dimensionless parameter was calculated and found to vary 
between 1.03 and 1.16. 

The seventh method used the Turbulent Boundary Layer (TBL) code and is shown Fig. 9. 
TBL method solves, simultaneously the integral momentum and energy equations for thin axi- 
symmetric boundary layers. 

The eighth method used backed-out ( h/C p ) data' that is used in CMA code to match the 
measured nozzle erosion profile. 

The following are observations and comparisons between the CFD results (Methods 1 and 2) 
and the other methods (Methods 3 to 8) as: 

(1) They all have similar profiles. They all increase to a maximum in the vicinity of the 
throat and then decrease along the nozzle wall with the exception of the correlation for the 
turbulent flow over a flat plate (Method 3). Method 3 shows the maximum to occur at the start of 
the nozzle (analogous to the leading edge of a flat plate) and decreases along the nozzle wall as 
expected. It was found that, along the last 40% of the nozzle length, the flat plate correlation 
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becomes reasonable. This is because curvature effects become less important along the nozzle. 

(2) Methods 1 and 2 (CFD), 7 (TBL), and 8 (backed-out) show the maximum to occur 
upstream of the nozzle throat. On the other hand, Methods 5 (Dittus-Boelter) and 6 ( Bartz ) show 
the maximum to occur at the nozzle throat. The latter group of methods is heavily dependent on 
the nozzle area ratio. For completeness. Method 4 (the modified Reynolds' analogy) also 
predicts the maximum to occur at the throat. 

(3) The profiles obtained from all Methods (with the exception of Method 3 (flat plate)), were 
analogous to the >> + distribution shown in Fig. 5. This confirms the dependency of heat transfer 
on the turbulence models as well as the value of v + . 

(4) At the nozzle throat: From the data used Fig. 9, the following values of ( h/C r ) are taken 
as 4, 4.32, 6, 12.09, 4.12, 3.89, 4.29, and 3.18 kg/m 2 -s for Methods 1 through 8, respectively. 
The values from Methods 1 through 8 differ Method 1 by 0.0%, 8.0%, 50.0%, 202.3%, 3.1%, - 
2.8%, 7.3% and -20.5%, respectively. This comparison is shown in Fig. 10. Results of Method 
4 fall out of the chart. Similarly, the values from Methods 1 through 7 differ from the backed-out 
data 23 (Method 8) by 25.8%, 35.8%, 88.7%, 280.2%, 29.6%, 22.3%, 34.9%, respectively. This 
comparison is also shown in Fig. 10. Results of Method 3 fall out of the chart. 

The comparison between Method 1 and Method 8 (backed-out) is within 26% and may seem to 
be in disagreement. The present results (Methods 1 and 2) are restricted to convective heat 
transfer. On the other hand, Method 8 (backed-out) was devised to match nozzle erosion that 
was a result of heat transfer (convective, radiative and conductive), chemical, and mechanical 
(particle impingement) effects combined. Therefore, comparing the present results and Method 8 
in a relative manner would be a sound approach. In the converging section of the nozzle 
(excluding the nose-tip), all methods over-predict the ratio h/C P . The heat transfer in the 
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convergence section may be inhibited by chemical effects. The nozzle wall is an ablating 
(transpiring) surface and was not modeled in this study. The boundary layer adjacent to the 
nozzle wall would be oxidizer-poor that inhibit chemical reactions and heat transfer. At the 
nozzle throat. Method 1 over-predicts Method 8 by 26%. It is to be noted that 12 data points 
used along the nozzle wall in Method 8. While, in the present study, there are 320 data used 
along the nozzle wall. All the methods over-predict the heat transfer has great implications. One 
concludes some factors that inhibit heat transfer. Some factors such as a thick boundary, layer 
and a transpiring surface would inhibit heat transfer. Furthermore, the condensed phase may 
flow as a film that would insulate the nozzle and inhibit heat transfer. In the diverging section of 
the nozzle, less agreement is obtained with Method 8 in the last 40% of the nozzle length. This 
disagreement may be attributed to eroded exit cone because of particle impingement based on 
observations made in post static and flight tests of this motor." Based on this author’s 
experience and observations, some conclusions have to be made through heuristic arguments 
alone. 

(5) In the vicinity downstream of the nozzle throat, the present results have a sudden 
decrease-increase-decrease. This sudden drop was attributed to the large drop in the specified 
surface temperature (Fig. 8). In the vicinity of the throat, the surface temperature dropped by 590 
K (1062 R) within 0.5 m (1.64 ft). This was easily verified by specifying uniform surface 
temperature at 2000 K (3600 R) and the sudden drop was not calculated. The heat transfer 
decreases when turbulent flow re-1 amin arize. Furthermore, this sudden drop was detected 
experimentally in the Jet Propulsion Laboratory ( JPL ) nozzle. 11, 1315 An isothermal wall 
boundary condition was used and was operated at a maximum chamber pressure of 1.72 MPa 
(250 psia). The JPL nozzle has a throat radius of about 2 cm (0.8 in) that amounts to about 3 
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percent of the throat size in this study. This is the reason for considering the possibility of 
turbulent boundary layer re-laminarization discussed in terms of acceleration parameter in the 
next section. 

Acceleration Parameter: The acceleration parameter, Ki, was calculated from" based on 

external boundary layer as 


K, = 


I r du ^ 


pM M. (*) 


V 


dx 


> 3x10 


-6 


(14a) 


and is shown in Fig. 11. The calculated values of K\ were smaller than the transition value of 
3x1 O' 6 . Therefore, re-laminarization of the turbulent boundary layer did not occur. The above 
acceleration parameter was translated by Coon and Perkins into terms more pertinent to tube 
flow as 




4 q w (jc) ( //(x) 

pJ(x)uJ(x)d(x)T{ C p 


> 1.5x10 


-6 


(14b) 


The maximum calculated values of K\ and K 2 were smaller than the transition values given above 
as 3xl0' 6 and 1.5xl0' 6 , respectively. Therefore, re-laminarization of the turbulent boundary layer 
did not occur. 


Summary and Conclusions 


The following summary and conclusions have been reached: 

• Two turbulence models and two types of wall treatment were used in this study. Based on 
economy and execution time, RNG k-£ with standard wall functions gives reasonable results in 
relative comparison with the backed-out data." 

• Convective heat transfer coefficients have been calculated using two methods. They have 
been compared with approximate methods, Bartz correlation, turbulent boundary layer (TBL) 
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theory code, and backed-out data based on post static and flight tests. 

• The interdependency of the convective heat transfer and wall y + has been shown. 
Therefore, consistency between a grid and a turbulence model cannot be over-emphasized. 

• The accuracy and the selection of the schemes and results are based on matching the RSRM 
ballistic predictions of mass flow rate, maximum head end pressure, and vacuum thrust and 
specific impulse and measured chamber pressure drop and was found to be good. 

• Convective heat transfer coefficients were calculated and compared favorably with existing 
theory. Qualitative comparison with backed-out data of the ratio of the convective heat transfer 
coefficient to the specific heat at constant pressure was made in a relative manner. Backed-out 
data was devised to match nozzle erosion that was a result of heat transfer (convective, radiative 
and conductive), chemical, and mechanical (shear particle impingement forces) effects 
combined. 
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Table 1: Summary of geometry, 

operating ballistic conditions, ballistic 
predictions, and gas thermophysical 
properties. 


Geometry 

r,/r* 

1.50 

rt/r* 

2.63 

Ao/A* 

2.065 

A,/A* 

7.72 

Ballistic Operating Conditions (TWR- 
16881. Rev. A fRef 18)) 

P „ (M Pa. tisia) 

6.28 (910.78) 

Chamber Pressure Drop (FSM-9 Static 
Test (Ref. 191) 

AP, (MPa, psia) 

1.138. 165 

Gas ThermoDhvsical Properties 

To (K, R) 

3419.2.6164.56 

fjo (kg/m-s, lbm/ft-s) 

9.25xl0’ 5 , 
6.22x1 0' 5 

C P (J/kg-K. Btu/lbm-R) 

1966.54. 0.47 

kn (W/m-K. Btu/hr-ft-R) 

0.397. 0.229 

MW (kg/kgmole) 

28.46. 


fl = jM, ( T/To / , k = k 0 (T/To) ; where £ = 


0.66687 as calculated from NASA-Lewis 
Program (Ref. 20). 



Table 2 CFD calculated pertinent results for the RSRM at 1 sec burn time 


Parameters 

Coarse Grid (105,850 cells) 

Fine Grid 
(376,100 
cells) 

Model 1 

Model 2a 

Model 2b 

Model 2c 

Viscous Model 

Standard k-f 

RNG a k-f 

RNGVf 


Near Wall Treatment 

Standard 

Wall 

Functions 

Standard 

Wall 

Functions 

Standard 

Wall 

Functions 

Two- 

Layer 

Zonal 

Model 

Desired wall y + 

V 

+ o 
o 

<N 

20 </< 
100 

20 < y + < 
100 

/< 1 

C,, 

0.09 

0.09 

0.055 

0.055-0.06 

Lit / 111 

10 4 

10 4 

10 6 

IliluiM 

I — — 

P S. max. CFD (M Pa, psia) 

6.29,911.90 

6.30, 913.62 

6.34, 

919.14 

6.28, 

911.16 

P o. max. CFD (M Pa, psia) 

6.29,911.99 

6.30,913.71 

6.45, 

934.92 

6.28, 

911.26 

% difference in P n (Using Table 1) 

0.13% 

0.32% 

2.65% 

0.05% 

M e. max 

3.55 

3.59 

3.58 

3.61 

% difference in mass flow rate 

-1.01% 

-0.92% 

-0.71% 

WmauMmi 

A P (M Pa. psia) 

1.38,200 

1.38.200 

1.23. 178 

1.17. 170 

% difference in A P (Using Table 1) 

21.21% 

21.21% 

7.88% 


Iterations used 

2000 

2000 

4000 

5000 


a :RNG (renormalization group theory), % difference = (calculated / predicted - 1) x 100 
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Fig. 2 RSRM cavity at prior to ignition time. 
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Fig. 1 Geometry and pressure distribution in the RSRM at 1 s post ignition time, 0, 0 < P s 
(M Pa, psia) < 6.34, 919.54; number of levels = 10; 0, 0 < APs (M Pa, psia) < 0.634, 91.95. 


29 




Model 1 (Table 2) 

Model 2a (Table 2) 

- - - * Model 2b (Table 2) 

— * -Model 2c (Table 2) 

■ QM-7 Test Motor [27] 

O QM-8 Test Motor [27] 

A FSM-9 [19] 






Fig. 4 Estimated height of cell adjacent to the 
wall 
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— r/r*, Nozzle radius — — • Model 1 (Table 2) 

- - - - Model 2a (Table 2) — - —Model 2b (Table 2) 


150 
125 
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75 V 
50 
25 
0 
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a: ir 



Fig. 5 Local wall jy + along the converging- 
diverging section of the RSRM nozzle at 1 s 
post ignition. 
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Fig. 6 Comparison of the velocity law of the 
wall with the RSRM nozzle calculations at 
exit. 
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Fig. 7 Comparison of the temperature law of 
the wall with the RSRM nozzle calculations 
at exit. 


35 







3 



Fig. 8 Local specified surface temperature and < 
post ignition. 





r/r 


r/r*, Nozzle Radius — — ■ CFD, Method 1 

- - - * CFD, Method 2 — - — Dittus-Boelter [17] - Method 5 

■■■ Bartz [10] - Method 6 ■ TBL Code [29] - Method 7 


O Backed-out [23] - Method 8 M Throat Location 



67.5 68.5 69.5 70.5 71.5 72.5 73.5 74.5 

x / r* 


Fig. 9 Local convective heat transfer on the nozzle wall of the RSRM at 1 s post ignition 
(Methods 3 and 4 are not shown). 
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Fig. 10 Comparison of heat transfer at the 


throat among all methods 
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Fig. 11 An acceleration Parameter ( K) based 
on externa] and internal boundary layers 
along the nozzle centerline of the RSRM at 1 
s post ignition. 
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